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Abstract. We derive an action-flux formula to compute the volumes of lobes quantifying trans- 
port between past- and future-invariant Lagrangian coherent structures of n-dimensional, transitory, 
globally Liouville flows. A transitory system is one that is nonautonomous only on a compact time 
interval. This method requires relatively little Lagrangian information about the codimension-one 
fSJ surfaces bounding the lobes, relying only on the generalized actions of loops on the lobe boundaries. 
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1. Transitory Systems. Finite-time transitions between steady states are com- 
mon in a wide range of physical systems. They play a key role in industrial mixing 
processes, mechanical systems in which parameters are modulated in a time-dependent 
^£^ manner, chemical reactions that progress to equilibrium, and shifts in local ecology 

^ due to a sudden environmental change. Since the transition mechanism may be com- 

Ch plex and the starting and ending states often differ, a prediction of the final state of 

i__i the system requires a detailed understanding of the transitional dynamics. 

In many cases, an analysis of transport and mixing in these systems can pro- 

K^ vide such an understanding. However, since any finite-time transition is aperiodically 

^_^ time-dependent, traditional techniques for computing dynamical transport p7 | [401 130 1 

^S| [26l[23] are often insufficient. In nonautonomous systems transport is often thought of 

OO as occurring between "Lagrangian coherent structures"; these are variously defined, 

CO for example, using ridges of finite-time Lyapunov exponent fields [15] HH [21], dis- 

Qf^ tinguished hyperbolic trajectories [T8[[19], or eigenfunctions of the Perron- Frobenius 

^^ operator [9l[T0]. However, few studies have quantitatively computed transport be- 

^^ tween coherent structures in aperiodic flows ^M, £, l32l [34], and these have been 

. . restricted to two-dimensions. Several studies of mixing in aperiodic flows have also 

^ been conducted [22l [38l [20] ; however, these have focused primarily on global mixing 

measures rather than transport between coherent structures, and again results have 

been restricted to 2D. To the best of our knowledge, no studies to date have given a 

C^ quantitative description of finite-time transport between isolated coherent structures 

in a 3D, aperiodic flow. 

In this paper, we present a formalism to compute transported volumes between 
Lagrangian coherent structures in a class of 3D aperiodic flows that we call "transitory" - 
the formal definition will be recalled from [34 in the remainder of this section. Our 
theory applies to incompressible vector fields that are, in addition, globally "Liou- 
ville", see ^ The Lagrangian coherent structures we consider are past- and future- 
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invariant regions of phase space, and transport between them corresponds to the vol- 
umes of certain lobes comprising the intersections of these regions. The salient point is 
that the computation of lobe volumes can be done by knowing only key "heteroclinic" 
trajectories that lie on the lobe boundaries. Compared to a naive, volume-integral 
approach, our method reduces by two the dimension of the Lagrangian information 
needed at any instant in time to compute a lobe volume. The result is an action-flux 
formula for n-dimensional lobe volumes, see ^ As examples, we will compute trans- 
port in a nonautonomous version of Arnold's ABC flow [3 in Q and in a model flow 
of a droplet in a microfluidic mixer in ^ 

On a phase space M, a transitory ODE [34] of transition time r has the form 

x = nx,t), nx,t) = {^W; m (1.1) 

where P \ M -> TM is the past vector field, F \ M -> TM is the future vector field, 
and ]/ : M X R ^ TM is otherwise arbitrary on the transition interval [0,r]. One 
way to model this type of behavior is by way of a transition function 

with the convex combination 

y(x, t) = (1 - s(t))P(x) + s(t)F(x). (1.3) 



For example, choosing 

s(t) = ^ (3-2-1 for ie[0,T], (1.4) 



t^ (. .t 



T 



along with (1.2), implies that V is C^ in time; we will use this form of s{t) in MJ 

Since the nonautonomous portion of the dynamics of ( |1.1[ ) is assumed to occur 
on a compact interval, it can be effected by a map. Suppose that V in ( |1.3| ) has a 
complete flow (^ti.to : M ^ M that maps a point from its position at t = to to its 
position at t = ti for any to,ti G M. Given a set Atg C M at time to, denote its 
evolution at time t under the flow by 

At = (pt,to{^to)i 

and its orbit in the extended phase space by 

A= {{At,t) :teR}CMxR. 

The orbit, A of any At^ C M is clearly invariant under (/p, and we refer to At as the 
time-t slice of A. The transition map T : M ^ M for ( |1.1[ ) is 

T(x)=(^,,o(x). (1.5) 

Consequently, a set Ao ditt = becomes Ar = ^(^o) at time r, and thereafter evolves 
under F. If the dynamics of P and F are known, then the only nontrivial work we 
must do is to characterize the map T. 

For transitory systems, it is natural to introduce some terminology for orbits 
according to their behavior under the stationary vector fields P or F. We will say an 
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orbit A is past invariant if At = Ao for alH < 0, and such a set is past hyperbolic if ^o 
is a hyperbolic invariant set of P. Similarly, A is future invariant if At = Ar for all 
t > r and is future hyperbolic if Ar is a hyperbolic invariant set of F. By extension, we 
call a slice At past /future invariant /hyperbolic if its orbit in the extended phase space 
satisfies the above definitions. Of course, sets that are invariant /hyperbolic under P 
or F need not be so under the transitory vector field V. For example, the orbit of a 
hyperbolic equilibrium p of P is both past invariant and past hyperbolic. However, 
because of the time-dependence of V on [0,r], T{p) is typically not an equilibrium 
of F, and even if it is, it need not be hyperbolic. Thus, the orbit of p under the 
transitory fiow (p need not be future invariant nor future hyperbolic. These concepts 
of "half-time" invariance and hyperbolicity will be used extensively in the remainder 
of the paper. 

Recah that stable and unstable sets ^'^'''(7) C M x R of an orbit 7 C M x R, 
are the sets of points that approach 7^ as t ^ +00 or —00, respectively. When 7 is 
past hyperbolic, each time-t slice of its unstable set for t < is the unstable manifold 
of the orbit of 70 under the stationary fiow of P; more importantly, the flow of this 
manifold under V is precisely the unstable manifold of the full orbit 7. However, the 
stable manifold of the orbit of 70 under the flow of P is not dynamically relevant for 
the transitory vector fleld — it almost certainly is not a stable set for 7. Thus, the 
unstable manifold of a past-hyperbolic set is dynamically relevant for the transitory 
vector fleld. Similarly it is the stable manifold of a future-hyperbolic set that is 
dynamically relevant for V. In the application described in Q the intersections 
between an unstable manifold of a past-hyperbolic set and a stable manifold of a 
future-hyperbolic set will be used to deflne lobes pivotal to the study of transport. 

2. Liouville vector fields. Recall that Hamiltonian systems are deflned on even 
dimensional manifolds M that are endowed with a closed, nondegenerate two-form 
uj G A^(M), the "symplectic form". A locally Hamiltonian [1, Prop. 3.3.6] vector fleld 
V : M X R ^ TM is one that preserves a;, that is, one for which 

jCv^ = 0, (2.1) 



where Cy is the Lie derivative (see (A.3) in Appendix [A|). Cartan's formula (A.4) and 



the assumption that duj = together give 

>Cya; = d{iv(^)' 



In this case, (2.1) implies that d{iYUj) = 0; in other words, whenever V is locally 



Hamiltonian, zycj is closed. 

If, in addition, this form is exact, 

zyoj = dH, (2.2) 

then V is globally Hamiltonian^ with the Hamiltonian function ilf : M x R ^ R. By 
Darboux's theorem [TJ Thm 3.2.2], there is a neighborhood of each point in M in 



which there are coordinates {q^p) so that uj = dq A dp. In these coordinates, (2.2) 
takes the form 

q = dpH, p = -dqH, 

i.e., a canonical Hamiltonian system. 



More generally, suppose an n-dimensional manifold M is endowed with a nonde- 
generate form, Vt G A^(M), i.e., a "volume-form". A vector field V is incompressible 
with respect to Vt^ or locally Liouville^ if 

Cv^ = (V • V)n = 0. (2.3) 

Liouville's theorem then implies that the volume of any region is preserved by the 
flow of V [31, §9.2]. As before, ( |2.3[ ) implies that iy^ is closed; if it is also exact, a 
global analog can be defined: 

Definition 2.1 (Globally Liouville [28^, §2]). A vector field V on a manifold M 
with volume form Q is globally Liouville if 

tv^ = d(3, (2.4) 

for some j3 G A^~^(M). Of course, if M has trivial cohomology then every closed 
form is exact, and there is no distinction between locally and globally Hamiltonian or 
Liouville vector fields. More generally, there may be some global obstruction to the 
existence of H or j3. For example, suppose that M = T^, and (^ = a; = dOi A d02 is 
the volume/symplectic form. Then the vector field V = pide^^ for a constant rotation 
vector p, is incompressible because iy^ = pid02 — p2d0i is closed. However, since this 
form is not exact (^i and O2 are not smooth functions on M), V is not Hamiltonian, 
or equivalent ly, not Liouville. 

As a second example, suppose M = R^ and Q is the standard volume, Q = 
dxi A dx2 A dxs- Using the natural identification of a two- form (" = CijkCidxj A dx^ 



with the vector ( = Qii and a one-form /3 = ^idxi with a vector /3 = /3^e^, (2.4) 
reduces to the statement that V = \/ x p. For instance, if F = Xidx^ is a Beltrami 
vector field on R^, i.e. V = \/ xV, then /3 = V, or 

/3 = Xidxi. (2.5) 

The ABC vector field (see Q is Beltrami, and hence has this property. 

The symplectic form u is, by definition, closed. If it is also exact, then there is a 
one- form u (often called the Liouville form) such that uj = —dv. Then if V is globally 
Hamiltonian, 

CyV = lydv + d{lyv) = d{lyV — H) = dL^ 

where L is the phase space Lagrangian. In canonical coordinates we could choose 
u = p ■ dq^ in which case L = p - q — H. 

Similarly, if a volume form Q is exact, we write 

n = da. (2.6) 

Then, if V is globally Liouville, 

Cya = lyda + d{iya) = d{zya + /3). 

Here zya + /3 G A^~^(M) is the Liouvillian analog of the Lagrangian L, so we make 
the following definition: 

Definition 2.2 (Lagrangian Form). Suppose that the volume form Q = da 
is exact and the vector field V is globally Liouville on a manifold M. Then the 
Lagrangian form A G K^~'^{M) is defined by 

Cya = dX , where A = tya -\- /3, (2.7) 
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and iv^ = dp. For example, the standard volume Q is exact on M = R^ with 

a = xsdxi A dx2. (2.^ 



When V = Xidj,. is Beltrami, (2.5), (2.7), and (2.8) imply 



A = esjkXsXj dxk + Xi dxi. 



(2.9) 



The discrete analogs of globally Liouville flows are exact volume-preserving maps; 
we define them here and note some key properties in anticipation of their use in ^ 

Definition 2.3 (Exact Volume-Preserving Map [24 ). Suppose that the volume 
form Q = da is exact on a manifold M. Then a diffeomorphism R : M ^ M is exact 
volume-preserving if there exists a generating form 77 G A^~^(M) such that 



a — R^a = drj. 



(2.10) 



Here we use the pushforward i?*, recall ( A.l ), instead of the pullback of [24 for later 
convenience. 

It is straightforward to show that if Ri and R2 are exact volume-preserving maps 
with generating forms 771 and 772, respectively, then the composition R = Ri o R2 is 
also exact volume-preserving [24 with generating form 



r] = r]i + i?i*7/2- 



(2.11) 



The generating form r] is the discrete analog of the Lagrangian form (2.7), and inter- 
acts with the latter as follows. 

Lemma 2.4. Suppose that Q = da is exact, V is globally Liouville with Lagrangian 
form Xv, and R is exact volume-preserving with generating form rj. Then the vector 
field W = R^V is globally Liouville with Lagrangian form 



A 



w 



R^X 



*Ay 



^WV- 



(2.12) 



Proof That W is globally Liouville follows directly from (A.5) and the invariance 
of Q under R^ : 



iw^ = R^iv^) = d{R^pv) := dp 



w, 



where ly^ = dpy Using (A.3) and Def. |2.2[ the Lagrangian form for W is derived 

by 



which gives (2.12). D 



As we will see in ^ the form A plays a central role in the computation of the 
volumes of lobes formed by the intersection of past- and future- invariant regions; it 
is analogous to the phase space Lagrangian we used to compute such volumes for the 
2D case [34 . 

3. Action-Flux Formulas for Lobe Volumes. In this section we will obtain 
the action-flux formulas to compute the transport fluxes. As a standing assumption, 
V will denote a transitory vector field ( |1.1| ) that is globally Liouville with respect to 
an exact volume form Q and that has a complete flow (pt,to- 
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Suppose that Vo, J^r ^ M are past- and future- invariant regions, respectively. 
By definition, trajectories within T^o at t = wih remain within it for ah t < 0: Vo is 
coherent under P in the Lagrangian sense. Similarly, Jv is coherent under the future 
vector field F. As a result, any transport between Vq and J> must occur during the 
transition interval [0,r], and the transported phase space itself is the collection of 
regions 1Zt = VtCiJ-'t, i-e., the intersection of V and T in any slice. We will call the 
components of the slices IZt, "lobes." Since V is Liouville, the intersection volume, 
or flux from Vo to Jv, 



^ = Yo\{VtnTt). 



(3.1) 
Vol(T(7^o)n 



is independent of time. In particular, using the transition map (1.5), ^ 

Let hi = dV and S = dT be the boundaries of the orbits of the future- and 
past- invariant regions so that the lobes, 1Zt^ are bounded by pieces of the slices Ut 
and St. As we will see below a key set in the action-fiux formulas will he X = U OS^ 
the set of orbits at the intersection of the lobe boundary components. 

In some cases U and S will be pieces of stable and unstable manifolds of future- 
and past-hyperbolic sets. For example, suppose that Vq and J> are topological balls 
whose boundaries are portions of the closures of the codimension-one unstable and 
stable manifolds of past- and future-hyperbolic equilibria, p and /, respectively!^ 
Under the transitory fiow, portions Ut C W^{p) and St C Wt{f) may intersect to 
bound a lobe 7^^, as sketched in Fig. |3.1[ a), which is also a ball. In this case, the 
intersection set. 



■UtHSt 



(3.2) 



is an (n — 2)-dimensional sphere. It is clear that W^{p) is past invariant and W^{f) 
is future invariant (though the lobe boundary surfaces Ut and St themselves are not), 
and, in this case, that Xt ^ po SiS t ^ — oo and Xt ^ fr sls t ^ oo. 





Fig. 3.1. (a) Simple lobe IZt C M formed by the intersection ofUt C W^(jp) and St C W^{f). 
(b) Similar to (a), except there are three surfaces that make up the lobe boundary: S^ , S^ C W^{f) 
andUt C W^(p). In both (a) and (b), intersections of the manifolds are shown as bold black curves 
and only the lobe boundaries dlZt are shaded. 



However, Xt need not be connected, p and / need not be single orbits, and IZt 
need not be a topological ball. For example, in Fig. 3.1 ^b), ft now represents a loop 



-•^For example, p could be the equilibrium p^ of the microdroplet flow of |5] see Fig. |5.2[ a), and 
thus Vo is the droplet itself. 



whose orbit is a future- hyperbolic periodic orbit, while p remains a past-hyperbolic 
equilibrium. In this case the lobe boundary contains two disjoint intersection curves, 
Xl and X|, both of which contract to po in the past, but in forward time approach the 
periodic orbit fr- As will become apparent, to compute the volume of 1Zt using the 
action-flux formulas, it is convenient (though not necessary) that the surface areas of 
Ut and St converge to zero in the appropriate limit in time. 

Since V is globally Liouville, the flux (3.1) is independent of time. It is this flux 



that we wish to compute, as it represents the portion of the past-invariant region Vo 



transported to the future- invariant region Jv. Stokes's theorem, using (2.6), allows for 



an immediate reduction of the volume integral for 1Zt to an integral over its boundary: 



<!> = Vol(7^^) 



Q 



7^t 



(3.3) 



dUt 



For the lobe IZt in Fig. |3.1[ a), dlZt =l^t^St^ while in Fig. 3.1 ^b) there are two disjoint 
surfaces S} ^ S^ C W^{f) on the boundary, so that dlZt =Ut -\- S} + S^. In general, 
dlZt can be decomposed into pieces Ul C dVt and SI C dFt of the boundaries of the 
regions Vt and J^t- Thus (3.3) becomes a sum of integrals of a over such submanifolds. 



It is to this computation that we now turn. 



Evaluation of (3.3) could be performed by numerical evaluation of the (n — 1)- 
dimensional surface integrals; however, this requires an accurate representation of the 
surfaces hit ^ St C dlZt^ which also implicitly requires knowing their time evolutions. 
This additional temporal information is essentially "wasted" since it is not explicitly 
used to compute the surface integrals. Furthermore, the exponential stretching typical 
of chaos can make obtaining well-resolved representations of hit and St computation- 
ally prohibitive. 

Our alternative reduces the dimension of the Lagrangian information necessary to 



evaluate (3.3) by computing the "generalized actions" of the orbits on the boundary 2"^, 



and requires evaluating an (n— 2)-dimensional spatial integral plus a temporal integral. 
In the extended phase space, this corresponds to an (n — l)-dimensional integral, but it 
explicitly uses the time evolutions of dhit and dSt that must be computed in any case. 
Our formulation applies to the case where dRt has components that are not stable 
or unstable manifolds of any future- or past-hyperbolic orbit (e.g., see ^. Indeed, 
since the orbit of any subset of M is invariant under (/:? in M x R, the result applies 
to general codimension-one submanifolds F^ C M. 

Theorem 3.1. Suppose that Vt = da is an exact volume form on M and Tt is a 
codimension-one slice of an invariant set of the flow cp of a globally Liouville vector 
field V. Then for any r gM, 



rt 



X\ ds -\- I a. 



(3.4) 



Proof. Differentiating a along the vector field V ^ employing Cartan's homotopy 



formula (A.4) (see Appendix [A|), and using (2.4) gives 

d 



dt 



a = Cyo^ = d{ivOi + /3) = d\. 



Integrating this expression from r to t using (A.5) then results in 



a 



^r,t<^ = / J~g^*s,t(^ds = / d{ipltX) 



ds, 



for any r. A second integration over Ft and rearrangement then gives 



/ a= / f / d{(pltX)jds^ ipl^^a 



dX \ ds -\- I a. 



which immediately reduces to (3.4) using Stokes's theorem. D 



It is interesting to note that the flow implicit in the time-integration in (3.4) 



may be chosen independently of the original transitory flow, provided it is globally 
Liouville. This observation will be used to simplify the computations in Qsj and an 
example of its implementation is discussed in Appendix [Bj 



Equation (3.4) simplifies if the surface area of Ft limits to zero in either backward 



or forward time. This typically occurs when Tt is a compact subset of an invariant 



manifold of a past- or future- hyperbolic orbit (e.g., as in Fig. 3.1). Under these 
assumptions, one can take the limit r -^ ±oo in Th. |3.1| to obtain the following. 

Corollary 3.2. Under the assumptions of Th. \3.1\ if the a-surface area ofTt 
vanishes as t ^ —oo, 



[ a= r ( [ X)ds, 



and if the a-surface area of Tt vanishes as t ^ +oo^ 

a = — I I / X] ds. 



dr. 



(3.5) 



(3.6) 



We will refer to (3.4)-(3.6) as the action-flux formulas. Since A is the n-dimensional 



analog of the Fagrangian, its integral along a codimension-two set of orbits gives a 
generalized action for that set. Thus, the action-flux formulas, in conjunction with 



(3.3), allow us to calculate the flux by computing the generalized action of sets of key 



orbits on the lobe boundary. 

For example, suppose that a lobe boundary d1Zr can be decomposed into N± 
connected submanifolds F^"^ that collapse in forward or backward time, respectively. 



Then (3.3), with (3.5) and (3.6), yields 



$ 



^7-00 vArr/'^^ hJr 



ar|- 



X]ds. 



(3.7) 



Implicit in (3.7) is a choice of orientation on the boundaries. We will always orient 



IZt with respect to a right-handed outward normal; this induces orientations on F-^ 
and FJ^, and, in turn, on their boundaries dTl~ and dT^^ . 

4. Example: Transitory ABC Flow. The ABC vector field [3], 



^Asinz + C cosy^ 

Bsinx -\- A cos z 

^C siny -\- B cosx j 



(4.1) 



models a steady, inviscid, incompressible Beltrami flow on T^. Interestingly, it is an 
exact solution to the Navier- Stokes equations, provided an appropriate forcing term 



is added to counter the effects of viscous dissipation; moreover for small Reynolds 
numbers, it is stable [11]. Despite the steadiness of the flow, its streamlines are 



chaotic [ini [Sol |T7] ; hence (4.1) is a prototypical example of a laminar vector field 
with complicated Lagrangian dynamics. 

The ABC vector field is locally Liouville; however, since there is no exact volume 
form on T^, it is not globally Liouville on this manifold. In order to apply the action- 
flux formulas of ^ we lift the z coordinate to R, letting the phase space become 
M = T^ X R. In this case, the standard volume Q = dx A dy A dz is exact on M, and 
we take a = zdx A dy. With this choice, (4.1) is Beltrami, and (2.5) and (2.9) give 



P = {A sin z -\- C cos y)dx + {B sin x -\- A cos z)dy + {C sin y 
X = {Asinz -\- C cos y) {dx + zdy) -\- {B sin x -\- A cos z) {dy - 
-\- {C siny -\- B cos x)dz. 



f B cos x)dz^ 
zdx) 



(4.2) 



Several recent studies have used finite-time Lyapunov exponents to analyze both 



steady and unsteady generalizations of (4.1) [131 [41]; these have focused primarily 
on extracting Lagrangian coherent structures by determining regions that experience 
maximal local stretching. Here, we study a transitory ABC flow in which the identi- 
fication of coherent structures of P and F is trivial, and focus on computing the flux 
between these structures. 



4.1. Transitory System. Modulating the coefficients A, 5, and C of (4.1) over 
a compact temporal interval makes it transitory in the sense of ( |1.1[ ). We choose to 
set C = for t < and 5 = for t > r, to give past and future vector fields 



Asinz 

P(x) = B sin X -\- A cos z 
B cosx 



F(x) = 



Asinz -\- C cos y 
Acosz 
C siny 



(4.3) 



The full transitory vector field is then given by the convex combination 

Asinz -\- s{t)C cosy 
V(x, t) = [ (1 - s{t))B sinx^A cos z 
s{t)C sin y -\- {1 — s{t))B cos x 



(4.4) 



as in (1.3), and we use the cubic transition function ( |1.4| ). We will denote the flows of 
P and F by cp^ and (/9^, respectively, and the flow of the full transitory vector field 



( |4.4D by (p. 

The autonomous vector fields P and F are integrable [8 ; indeed, they have in- 



variants 



Hp (x) = Bsinx -\- A cos z, 
Hf (x) = A sin z -\- C cos ?/, 



(4.5) 



respectively. Moreover, these functions act as Hamiltonians that generate the flows of 
P in the (x, z) plane and F in the (?/, z) plane. Since y = Hp for P and x = Hp for F, 
the motion in these transverse directions is trivial; consequently, the flows of P and F 
can be completely characterized by their two-dimensional portraits, see Fig.|4.1[ Note 



that the level sets of the invariants (4.5) become two-tori in M, with the exception of 



certain critical sets, which correspond to separatrices. 
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Fig. 4.1. Poincare sections of the invariant two-tori of (a) P and (b) F for (A,B,C) = 
(1,0.3,1.5). The invariant manifolds of the equilibrium p of P are the bold (red) curves in (a), and 
those of the equilibria f^ of F are the bold (blue) curves in (b). The ranges of {x,y,z) are shifted 
to visualize the resonance zones Vo and J^^ . 



In each 27r range of z, P and F generally have two elliptic and two hyperbolic 
periodic orbits; in the special case A = B = C = 1^ these become lines of fixed points. 
We will set A = 1, without loss of generality, and assume that 



0<B<A = 1<C; 



(4.6) 



with this choice, the resulting phase portraits are like those in Fig. 4.1 The past- 
and future-invariant regions that we will analyze are bounded by the manifolds of the 
hyperbolic periodic orbit 

p={(f,y,7r) I ye[0,2Tr]} 

of P and the hyperbolic periodic orbits 

f = {{x, TT, f + 27r(fc -I)) \x€ [0, 27r]}, (4.7) 



of F (two of these are shown in Fig. |4.1[ b)). 

The invariant manifolds — under P — of p correspond to the level set 



Hp{^) 



H^P 



B-A. 



(4.^ 



Because of the horizontal periodicity of M, W^ip) forms a pair of homoclinic connec- 
tions, shown in cross-section in Fig. |4.1[ a). Each is homeomorphic to a two-torus in 
M, as shown in Fig. |4.2[ and together, these manifolds bound a past- invariant region 
Vq. Similarly, the invariant manifolds — under F — of /^ correspond to the level set 



Hf{-x) =H^:=A-C. 



(4.9) 



Sinc e M is unbounded in z, these form heteroclinic connections, see Fig. |4.1[ b) and 
Fig. 4.2 'a). The heteroclinic connections between /^ and /^^^ bound a future- 
invariant region Tl^; for each k e Z these are just shifted copies of F^, shown in 



Fig. 4.2 



It is important to remember that Vo and T!^ are not invariant under the transitory 
flow cp. However, they are Lagrangian coherent structures, or "resonance zones," of 
the past and future vector fields, respectively. We consider the problem of computing 
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Fig. 4.2. (a) The past- and future-invariant regions for \4--4\ with (A,B,C) = (1.0,0.3,1.5) 
are bounded by the invariant manifolds W^(p) (red) and W^(p~^) (blue), respectively. The lobe 
shown corresponds to IZr = Vt H J-"^ for r = 0. Its boundary consists of three surfaces lAr C W^{p), 
S]. C W^{f'^), and S^ C W^{f-^), which intersect in two loops X^ (green) shown in (b). 



the transport from Vq to each of the J^^; that is, the volume of the lobes IZ^ = Vr^F^. 
There is always at least one such lobe for any r and choice of parameters subject to 
( |4.6[ ), and when r is finite there are only finitely many. The accompanying movie file 
(ABC Flow Lobes Vs. Transition Time) shows the lobes at t = r, for increasing 



values of r and for the same parameters as in Fig. 4.2 



4.2. Computation. The volumes of the lobes IZ^ will be computed using the 
action-fiux formulas (3.5) and (3.6), which rely on knowing the orbits of the intersec- 
tion curves 



x, = M/;(p)n|Jl^^(/'=). 



kel 



As an example, two such curves, Z^ and Z^, are shown in Fig. 4.2 'b). For the moder- 
ate values of r that we study below, only the intersections of Vr with J^l and J^^ are 
nonempty, and thus the only lobes formed are IZ]. (the primary lobe) and 7^^ (the sec- 
ondary lobe). To simplify notation, we adopt the convention of referring to elements 
of the secondary lobe with a "tilde" (i.e. 7^^-, X^-, etc.), and omit the superscripts for 



both lobes (cf. Fig. [45). 

The computation of the intersection curves is done with a root finding and con- 
tinuation method, and is simplified by using the level sets (4.8) and ( |4.9| ). It is 
convenient to parameterize the past manifold as G : T^ ^ Wq{p) = dVo, and 
search for intersections in parameter space. Using (u^v) as the parameters, (4.8) 
gives G{u^v) = {x{v)^u^ z{v))^ where 



x{v) 
z{v) 



2v 



ve [0,7r), 



^-2v, vG[7r,27r], 



cos 



-i(2fsin2t;-l), 



27r-cos-i (2f sin^u-l) 



v€ [O.tt), 
V e [tt, 27r], 



(4.10) 



see the sketch in Fig. 4.3 Note that increasing v corresponds to a counterclockwise 
circuit around the separatrix loop in the xz-plane, while u is simply the ^/-coordinate. 
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T2[0,27r] 

27r| ; ,p 



^ 



G(u,v) 




Fig. 4.3. Parameterization of Wq(p). Arrows along the manifold are included to indicate 
orientation. 

The parametric representation of the intersection curves is then given in terms of 



the level set (4.9) by 



lG = {{u,v)eT^ I Hf{T{G{u,v))) = H*f}, 



where T is the transition map (|1.5|). Two examples are shown in Fig. 4.4 



see 



Appendix [B] for a discussion of the continuation techniques used in the computa- 
tion of these curves. The corresponding intersections in phase space then become 
^t = ^^,o(G(Xg')), and the curves shown in Fig. 4.2 'b) are the phase-space represen- 
tation of those in Fig. 4.4 'a) under this mapping, with t = r = 0. 




Fig. 4.4. Intersection curves Xq in parameter space for (A, B, C) = (1, 0.3, 1.5) with r = (a) 
and T = 2.0 (h). The accompanying movie file (ABC Flow Manifold Intersections Vs. Transition 
Time) shows the dependence of Xq on transition time, for r G [0.5,7.5]. 

Given numerically-computed curves Xg, it is straightforward to compute the flux 
$ using the action-flux formulas. We first discuss the case where there is one lobe, 
the prima ry lo be IZr = Vr^J^l^ and its boundary has the form dlZr =V(rUS^US^ 



as in Fig. 
given by 



4.2 



'a). This occurs when B < ^ and r is small enough. The flux is then 



^ = Vol(7^^) =/ a+/ a+/ a. 



(4.11) 
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Note that under the flow (/:?, S^ -^ /^ and S^ -^ f^ as t ^ +00, and moreover, since 
these surfaces he on the local manifolds W^{f^)^ their surface areas limit to zero. 



Thus, by (3.6) 



7<S:^ Jt \Jxl 



X ds. 



(4.12) 



where X^ = dS^^ and the orientations of these contours are chosen consistent with 
a right-handed outward normal to the lobe. Further note that over the range of 



integration in (4.12), Lp = Lp . It is also possible, and may be more efficient, to 
compute these integrals backward in time under Lp^ (see (B.2) in Appendix [b|) . 



To compute the integral over Ur in (4.11 ) we first use ( 3.4[ ) to pull the integration 
back to t = 0: 



a 



X\ ds -\- I a. 



(4.13) 



The ffist integral on the right-hand side is easily computed numerically, as it is over the 
compact transition interval. Computation of the second is a bit more difficult: since 
Uq encircles Wq{p)^ it does not collapse top under (/9^, and consequently, its a-surface 
area does not vanish in either direction of time. To get around this, we can divide Uq 
into subsurfaces that do collapse under (/9^, reducing the last integral in (4.13) to a 



sum of integrals over these subsurfaces, each of which can be evaluated using (3.5)- 
( |3.6[ ) (see (B.l) and technical remarks on this splitting in Appendix [B|) . With these 
techniques, the flux (4.11 ) can be computed using (4.12) and (4.13). Several additional 
techniques can be used to speed up the computations and decrease numerical errors, 
see Appendix [B] 

If 5 > ^ , the past-invariant set Vr intersects both J^l and J^^ to form a primary 
lobe IZr and a secondary lobe IZr^ even when r = 0, see Fig. 4.5 A secondary lobe 




Fig. 4.5. Lobes at t = r for (A, B, C, r) = (1, 0.8, 1.5, 0). The dotted green segment is added to 
dSr and dS^ to ensure that these are closed curves. 



also exists when B < ^ 



4.4 



provided r is large enough; for example, for the parameters 
IZr is formed at r 



4.5. This lobe can be seen in Fig. 



4.6 



for 



used in Fig. 

r = 6. The total flux ^ is then simply the sum of the volumes of the two lobes, and 
computing each is similar to the single- lobe case described above. 
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Fig. 4.6. Intersection of Vt with T^ for gTTp with {A,B,C) = (1.0,0.3,1.5) and r = 6. 
Intersection curves are shown (a) with their corresponding invariant manifolds W^(p) (red) and 
H^^(/0'i'2) (blue), and (b) without them. 



4.3. Results. Using the techniques discussed above, we computed $ as a func- 
tion of the parameters B and r. The volumes of the primary and secondary lobes, 
each of which contributes to <l>, were computed independently. A summary of these 
lobe volumes as percentages of the volume of the past- invariant set Vq is shown in 
Fig. 4/7 The volume of Vo itself is given by the quadraturq^ 



Vol(7^o) = 47r 



2^^ 



cos ^ (§(1 — sinx) — l) dx 



(4.14) 



The curves in the figure denote the parameters at which the second lobe IZr emerges. 
Thus for B values below this curve, the volume of the secondary lobe in Fig. |4.7[ b) is 
zero. Note that for the parameter ranges of the figure, the fiux due to the secondary 
lobe is never more than 3.5% of the volume of Vo- Also, the percent fiux of the 



primary lobe does not strongly depend upon r; this variation is shown in Fig. 4.8 for 
B = 0.5. 





Fig. 4.7. Ratio of the volumes of the (a) primary and (b) secondary lobes to Vol(7^o)- Here 
(A,C) = (1.0, 1.5)j and B and r vary. The curves denote emergence of the secondary lobe. 



^Similarly, the volume of >F^ is given by the same formula upon replacing ^ with -^ and x with 
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28 

?27.5 



^O 26.5 

> 



- Sum 
■ Primary 



Fig. 4.8. Dependence of the primary and secondary lobe volumes on transition time r for 
(A, B, C) = (1, 0.5, 1.5). The dotted line denotes the primary lobe volume and the solid line denotes 
the total flux <l>. 



We performed two checks on the accuracy of the numerical integrations. For 
: 0, the intersection curves Xq = ^r correspond to the points 

z = cos~"^ (z(-'- ~ sinx) — l), 



y = cos ^ (^(1 — sinz) — l), 

as X ranges over [f , ^]. Since analytical expressions exist for both the intersection 
curves and the lobe boundary itself, it is straightforward to numerically compute the 



lobe volume directly from the 2D integral in (3.3). These computations revealed that 



the relative error in the action-flux computations was, on average, on the order of 

As a second check, and for nonzero r, we estimated ^ using a Monte Carlo 
simulation. We uniformly seeded Vo with N points, advected each point to t = r, 
and determined the number Ni^ of orbits of sample points that lay within the J^^ at 
t = r. An estimate for ^ is then 



$ 



MC 



Vol(7^o) 



N 



using (4.14). The relative error in any realization of this Monte Carlo computation is 



estimated by [Ml §7.6] 



1 / N 



1. 



(4.15) 



For N = 10^, the difference between the action- flux computations and the Monte 



Carlo simulations was typically less than (4.15). For larger values of r the difference 



between the two methods increases; this is likely due to resolution issues with the 
intersection curves X,-, as neighboring points along these curves begin to separate 
significantly as r becomes large. 

Note that a Monte Carlo computation of the flux is only feasible due to the 
simple nature of Vq and J^^. Since the boundaries of these past- and future- invariant 
sets are known analytically, it is straightforward to both uniformly sample Vq and 
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determine which trajectories he within F^ dX t = r. This could be computationahy 
prohibitive for more comphcated past and future- invariant sets. In such cases, using 
the action-flux formulas to compute ^ typically requires less Lagrangian information. 

5. Example: Microdroplet Flow. As a second application of our theory, in 
this section we study transport between two halves of a droplet moving through a ser- 
pentine microfluidic channel mixer. Microfluidic devices have numerous applications, 
for example to detect specific antigens in the blood [7 , perform macromolecular or 
cellular assays and analyze DNA [4 , and even filter circulating tumor cells from the 
blood for early stage cancer diagnosis and metastasis detection [35j|43]. They have 
also been key components of process intensification in industry [3T, ^29^ , an effort to 
decrease processing time, make more efficient use of raw materials, and gain greater 
functionality and sensitivity with respect to device size. 

Many applications of microffuidics require thorough mixing; however, small length 
scales or high ffuid viscosities often force a Stokes ffow regime in which mixing by 
molecular diffusion alone can be impractically slow [36 . Consequently, chaotic advec- 
tion is required [2^, TO , and designing devices in which this occurs is of much interest 

We model a channel mixer for which the mixing occurs within a droplet formed 
by the injection of equal volumes of ffuids A and 5, see Fig. |5.1[ a). For simplicity we 
assume that the droplet is a sphere and that each ffuid occupies one hemisphere at 
t = 0. The plane separating A and 5 at t = is the injection plane, denoted by Uq. 
Subsequent motion of the droplet through the serpentine microchannel, as sketched 



in Fig. 5.1 'b), will induce a time-dependent vector ffeld ^(x, t) within the droplet, 
with the goal of stretching and folding the initial interface to enhance the mixing by 
molecular diffusion. Finally, after a transition time r, the droplet is extracted by 
dividing it into two hemispheres along an extraction plane, Sr- 



(a) 

X 

A 



-^z 



Uo = {x = 0} 





t = 



Fig. 5.1. (a) Streamlines in the xz -plane (comoving frame) and initial locations of the fluids A 
and B and the injection plane Uq that separates them, (h) Streamline dependence on channel shape. 



We assume Uq is invariant under P(x) = F(x, 0) so that the hemispheres A and 
B are each past invariant. Similarly, we assume Sr is invariant under F(x) = V(x,r). 
Extending the dynamics to t G M in this way gives a ffow that is transitory in the 
sense of ( |1.1| ). Our goal is to compute the fraction of A in each of the extracted 
hemispheres as a function of the microchannel shape. 

For simplicity we assume that the droplet remains spherical as it moves through 
the serpentine channel. We also assume that it is always in contact with the channel 
walls, so that — in the lab frame — the velocity at these contact points must be zero. 
For a straight channel, the resulting steady ffow is axisymmetric with a vortical recir- 
culation, such as that sketched in Fig. |5.2[ a), and this is conffrmed experimentally [42^ 
for a "bullet-shaped" droplet in a straight channel with rectangular cross-section. It 
is reasonable to think that these results would extend to a droplet that has sufficient 
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surface tension to maintain a spherical shape, as weh as to a more symmetric, circular 
channel. Finally we assume that the center of mass velocity of the droplet remains 
parallel to the channel walls so that, as the channel bends, the droplet's velocity field 
rotates to keep its axis of symmetry parallel to the walls, as sketched in Figs. |5.1[ b) 
for 2D, andlS^b) for 3D. 




Fig. 5.2. (a) Streamlines for the droplet vector field Vq in the xz- (blue) and yz- (red) plan es, 
(h) Corresponding contours of the rotated vector field with angles 0, ip and cf) as defined by ^5.5}^ . 

5.1. Stationary Model. The steady, low- Reynolds-number fiow inside a spher- 
ical droplet subject to a uniform external creeping fiow \J = —Uz (the Hadamard- 
Rybczynski problem) is known analytically [12j. Following ^46], we take this fiow 
to model the steady motion within a droplet moving through a straight section of 
microchannel. To nondimensionalize, we normalize lengths by the droplet radius, D^ 
velocities by /7, and time by 4D(1 + l~j)/U ^ where ji is the ratio of the viscosity of the 
fiuid within the droplet to that of the exterior fiuid. The nondimensional vector field 
within the droplet in a comoving reference frame is given by 



Vo(x) 



,2(1 



2xz 

2yz 

2x^ - 2y2 



^'). 



In this frame, the droplet sits at the origin and its boundary (r^ = |x|' 



1) is 



invariant under Vq. Certainly the length scales in typical microfiuidic devices imply 
that the low Reynolds number assumption used to obtain this solution is appropriate. 
In addition, if the surface tension at the droplet boundary is sufiiciently large, the 
assumption that the droplet remains spherical seems reasonable. 

The vector field Vq has two hyperbolic equilibria (stagnation points) p^ = (0, 0, 1) 
and p^ = (0, 0, —1), both on the droplet boundary, recall Fig. 5.2 ^a). The correspond- 
ing two-dimensional manifolds W^{p^) and W^{p'^) form a heteroclinic connection be- 
tween p^ and p^ that comprises the droplet boundary. Similarly, the one-dimensional 
manifolds W^{p^) and W^{p'^) form a second heteroclinic connection between these 
hyperbolic points in the interior along the axis of symmetry. There is also a ring of 
elliptic equilibria at z = with x^ + 7/^ = 1/2. 

The vector field Vq can be viewed as a sum of two, 2D Hamiltonian vector fields 
with (x, z) and (y, z) as canonical variables, and Hamiltonians 



ff(x) = 


= (1- 


-x'- 


-y'- 


- z^)x, 


K{^)-- 


= (1- 


2 

- X - 
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-y'- 


-z^)y, 



(5.1) 



respectively. Then Vq is equivalent to 

Vo(x) = 




(5.2) 



With this formulation, it is straightforward to show that (5.2) is globahy Liouvihe in 
the sense of Def. 12. 11 with 



Po = H dy — K dx^ 

Ao = {zx -\- H) dy — {zy + K) dx. 



(5.3) 



The vector field Vq contains no swirl: each plane containing the 2:-axis is invariant. 
Axisymmetry implies that the dynamics in each such plane is equivalent to that in 
the y = plane, for example, which is Hamiltonian with H{x, 0, z). Consequently the 
dynamics of Vq is completely integrable. 

5.2. Transitory System. In our model, as the droplet moves in the serpentine 
channel, e.g.. Fig. |5.3[ transitory time dependence is introduced when the vector field 



(5.2) rotates to maintain its axis of symmetry parallel to the channel walls. 



(a) 



MicroChannel Wireframe 



Centerline: Side View 




Fig. 5.3. (a) Wireframe of the microchannel for angles given by \5. 1 0\ with r = 3.5 and 
^ = 7r/8. The red centerline is equidistant from the four corners of each cross-section and everywhere 
tangent to the bulk velocity of the droplet, (b)-(c) Projections of the red centerline onto the xz- and 
yz-planes. 

We will assume that the fiuid A initially occupies one of the hemispheres bounded 
by the injection plane Uq- Without loss of generality, we will choose Uq = {(0, ^^ ^)}|j 
and suppose that A occupies the "negative" hemisphere (x < 0), denoted by Vq- Note 
that Vq is past invariant under V, since Vt = Vq for t < by assumption. Our goal 
is to compute the fiux ^ of A from Vq to one of the hemispheres bounded by the 
extraction plane Sr- We will investigate two possible choices. 



Sr = {(0,y,z)} or Sr = {(x,0,2:)}. 



(5.4) 



•^Note that here Uq is not the unstable manifold of any past-hyperbolic orbit; indeed, its surface 
area does not vanish in either direction of time. We nevertheless use a notation consistent with that 
of the previous sections, i.e. Uq is a past invariant surface rather than the past unstable manifold. 
Similar notation is used for the future-invariant surface Sr- 
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Let Tr denote the "positive" extraction hemisphere; i.e., x > or ?/ > 0, for the cases 
above, respectively. Note that the hemisphere J> is future invariant. Since volume 
is preserved, the flux we compute is equal to that of the fluid B to the complement 
of Tr- The values of the fluxes of A to the complement of J> and B to J> follow by 
volume preservation. 



The flux, of course, depends heavily on the choice (5.4) of extraction plane. For 
example, if r = and S^ = %, Tr will not contain any of fluid A\ <l> = 0. However, if 
the angle between Sr and lA^) is |^, as for our second choice, A and B will be equally 
represented in each extracted hemisphere when r = 0, so <l> = 7r/3, one-quarter of the 
droplet volume. 

To model the shape of the serpentine channel, we will specify the rotations that 
give its axis and orientation at each time t. Let d(i) and '^(t) be the angles between 



the channel axis and its projection onto the lab-fixed xz- and ?/z-planes, see Fig. 5.2 



The angle ^(t) will represent an additional torsion about the axis of symmetry. These 
"Tait-Bryan angles" correspond to the matrices 



Ry{Q) = 



cosO 








1 


-sin^ 






sin^ 


cos^ 



,R.W = 



1 










COS(f) 


— sin (/) 





COS'0 


sin-^ 


,i?.(0) = 


sin(/) 


cos(/) 





— sinV^ 


cosV^ 







1 












(5.5) 



and give an overall rotation 



R{t) = Ry{e{t))R,{iP{t))R,{<P{t)). 



(5.6) 



These are not the typical "right-handed" rotations in R^, but are designed to give 
the horizontal, vertical, and torsional angles along the microchannel. In the frame of 
reference of the center of mass of the droplet, the time-dependent vector field is then 



the pushforward, (A.2), of Vb by R{t)^ 



F(x,t) = i?*(t)Fo(x). 



(5.7) 



Since the channel is initially aligned with the z-axis, 6>(0) = '0(0) = 0(0) = , and 
since the process only occurs for t G [0,r], we can formally extend the vector field to 
t G M by setting 

e(t) = i^(t)= (j){t) = for t^ [0, r] , 



so that (5.7) is transitory. In similar fashion to ^5.1, denote the hyperbolic equilibria 



at the leading end of the droplet by p^ and /^ and those at the trailing end of the 
droplet by p'^ and /^, under the past and future vector fields, respectively. Since the 
droplet boundary at r = 1 remains invariant under (5.7), the phase space M is simply 
the closed ball of radius 1, centered at the origin. 

Since the droplet is assumed to contact the channel walls, no-slip boundary condi- 
tions and ( |5.2[ ) imply that, in the lab frame, the droplet moves with a nondimensional 
speed of two. Since its velocity is aligned with the channel axis, the position of the 
droplet center, c(t), satisfies the initial value problem 



c{t) = 2R(t)z , c(0) = 0. 



(5.^ 



The function c(t) also defines the channel center as a function of the nondimensional 
time t; moreover, since the dimensional arclength along the channel is simply s = 
AD(1 + /i)t, c{t) also defines the channel center as a space curve, recall Fig. 5.3 
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Note that from (5.6), the torsion (j) drops out of (5.8), as it should. Torsion, 
however, does affect the channel shape. Indeed, assuming that the channel maintains 
its shape in the direction perpendicular to the current channel axis, c(t), then a point 
w(0) on the channel wall at the injection point evolves to w(t) = c(t) + i?(t)w(0). 
In order that this model correspond to a physically realizable channel, the walls must 
correspond to an embedded submanifold — there can be no self-intersections. For a 
circular cross-section, self-intersections will occur if the local radius of curvature of 



any bend in the channel is less than the cross-sectional radius. For rotations (5.6), 
a circular channel with a nondimensional radius of one will be physically realizable 
when 

i92(t)cos2V^(t)+^2(t) <4. (5.9) 

For example, it is sufficient that both |^| and \ip\ < 2. If the cross-section is not 
circular, it is more difficult to obtain an analytical requirement. 
For the computations below, we will use the shape functions 

e{t) = ilj{t) = ^sin(27rt/r) , ^(t) = 0, (5.10) 



for t G [0,r], recall Fig. |5.3[ and fix 0{t) = ilj{t) = (j){t) = outside the transition 
interval. Since these angles vanish at t = and t = r, both P and F are equal to 
Vq in ( |5.2[ ). For this model, there is a critical transition time r* for each amplitude 
^, below which the channel walls would self-intersect. For a circular channel this is 



easily obtained from (5.9): 



When the channel has a square cross-section, as shown in Fig. |5.3[ a slightly larger 
r* is required because the self-intersection first occurs at the corners; a numerical 
solution for two cases gives 

r 2.4675, e = f 
'^ \ 4.9348, ^ = f ' 

In terms of dimensional variables, the physical channel length is L = [/r, and so, for 
a given viscosity ratio, these requirements imply a lower bound on on the "inverse 
aspect ratio" of the channel, 

^>4(1+m)t*. 



Even when r < r*, the model (5.7) still corresponds to a droplet that is forced to 
rotate; however, to realize these rotations in a lab, some other mechanism, such as 
electromagnetic manipulation of a charged microdroplet [25l |45] would have to be 
used. 



We will use the model (5.10) to investigate the effects of the transition time r and 
amplitude ^ on the transported fiux. As we noted above, given a fixed diameter and 
viscosity ratio, the transition time r is a proxy for the channel length. The oscillation 
amplitude ^ controls the magnitude of the bends in the channel. 

5.3. Computation. To compute <l>, we employ the action-fiux results of ^ For 
any r > 0, there exist well-defined lobes containing all the fiuid A in Jv; these lobes 
are bounded at any time by slices of the orbits of the injection and extraction planes, 
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hit and 5^, and the invariant droplet boundary dM (see Fig. 5.4). Since Uq and Sr are 
not hyperbolic manifolds, the intersections of their orbits in any time-t slice are not 
heteroclinic, as was the case in Q however these intersection curves are still needed 
for the computation of the lobe volumes according to the action-flux formulas. They 
are interior to the droplet (except at possibly two points on dM) and we denote them 
at time r by 

x^* = T(%)n5,. 

where T is again the transition map \l.b\ . The intersections of T{Uq) and Sr with 
dM itself are also important. They lie on the invariant spherical boundary and are 
denoted by 

X^ = T{dUQ)\JdSr. 



T = l 




r = 2 




Fig. 5.4. Lobes contributing to the flux ^ at t = r for Sr = {x = 0} and ^ = 7r/4. The 
intersection curves X are shown in green and the spherical boundary dM is not shown. 

We use a method similar to that in Q to identify the intersection curves. The 
injection plane Uq has a natural parameterization G : ^i(O) -^ M, given by 

G{u^v) = {O^u^v). 

Letting W : M ^ R so that the extraction plane Sr is its zero level set, the interior 
intersection curves in parameter space. 



T-int 



{{u,v)eV,{0) I W{T{G{u,v))) = 0}, 



(5.11) 



are computed using a 2D root finder. Figure |5.5| shows these curves for ^ = 7r/4, 
various values of r, and both choices for the extraction plane (5.4). The corresponding 
curves in phase space at t = r are 

X;"*=T(G(X^"*)). 

To compute the transition map, we use a combination of Cartesian and spher- 



ical polar representations of (5.7) and integrate the vector field using MATLAB's 



implementation of the Runge-Kutta (5,4) Dormand-Prince pair. The point is that in- 
tegration in Cartesian coordinates does not respect the invariance of the sphere 9M, 
while standard spherical coordinates induce singularities at the origin and along the 
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[5. 1 jp with ^ = 7r/4 for St = {x = 0} (top row), and Sr 
p to the 



Fig. 5.5. Intersection curves Xq^ 
{y = 0} (bottom row). Regions that map to the positive side of the extraction plane (poi nts for which 
sgn(yV(T(G(w,'u)))) > 0^ are shaded; these correspond to the subsurfaces ofUr C M in (5. j^p along 
whose boundary curves we must integrate to compute <l>. Representative orientations are shown in 
the left column. The accompanying movies Micro drop let Intersection Curves: Extraction Plane 
{x=0} and Micro drop let Intersection Curves: Extraction Plane {y=0}, corresponding to the top 
and bottom rows, respectively, show the continuous dependence of Xq^ on r. 



positive 2;-axis. To avoid the latter, we introduce a second spherical representation 
in which the inclination is measured from the positive ?/-axis. We monitor both incli- 
nation and radius, switching between spherical representations when the inclination 
drops below some prescribed minimum, and switching to the Cartesian representation 
near the origin. 

We turn now to the application of the action- flux formulas to compute ^. It is 



generators, (2.10), 



not hard to show that the rotations in (5.5) are each exact volume-preserving, with 



r]y = -cos sin 0{z —x)dy — ysm 0{z dx -\- xdz)^ 
Vx = - cos'0sin'0(y^ — z'^) dx -\- xsin^ ip{zdy + y dz)., 



respectively. Consequently, R{t), (5.6), is also exact volume-preserving and its gener- 
ating form, by (2.11), is 



v{t) = Vy{t)^Ry4t)vAt)' (5.12) 

Finally, since Vo is globally Liouville, Lem.|2. 4] implies that the transitory vector field 



V(t), (5.7), is as well, with the Lagrangian form 



X{t) = R4t)Xo^Cvit)V{t). 



as obtained from ( |2.12[ ) using ([53|, ( [5^ , ( [5J| ), and ( |5.12| ). 

To apply the action-fiux formulas, we must identify the appropriate boundaries 
over which to integrate. There are two possible types of lobes (interior lobes 1Z\^^ 



and boundary lobes T^f, see Fig. |5.4[) and three possible types of bounding surfaces. 



corresponding to portions of Z//, S and dM. We need only specialize (3.4) to integrate 
the 2-form a over these bounding surfaces. 
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Since the integral of the two- form a over Uq is identicahy zero, (3.4) imphes that 
an integral over Ur can be reduced to 



m 



dUi 



A ds. 



(5.13) 



A similar simplification applies to Sr^ for the two cases that we study. 

The integrals over portions of dM can be simplified by noting that dynamics on 
this surface under the stationary vector fields P = F is trivial due to axisymmetry 
and the invariance of dM. Indeed, if p and ^ are the azimuthal and inclination angles, 
then the fiow of F reduces to 



(pf(p,C) = (p,2arctan(tan(C/2)e2(*-*°))) 



(5.14) 



on the droplet boundary. We use this result to compute the action integrals outside 
the transition interval. 

Since dM is heteroclinic from pi to /2 under the transitory fiow, an integral over 
a surface Br C dM can be computed using the fiow of F in either (3.5) or (3.6): 



Jb^ 




A ds 



Ad^^) 




A ds. 



Ad^A 



(5.15) 



Here we simply chose to evaluate the integral that converges faster. It may happen 
that one the equilibria /* of F lies on the boundary dBr] in this case, only one of the 
two integrals (5.15) converges. In the rare case that both equilibria lie on dBr^ we 



can modify the fiow F by applying a rotation so as to effectively move the equilibria. 
Convergence of these integrals can also be accelerated by estimating their exponential 
tails as discussed in Appendix [B] 

As always, care must be taken to ensure that the intersection curves are oriented 
to be consistent with a right-handed outward normal to the lobes. For example, the 



correct orientation for integration over dlA is indicated in Fig. 5.5 



5.4. Results. Figure 5.6 summarizes the fiuxes computed using the shape model 



(5.10) for two values of the amplitude ^, the two choices of extraction plane in (5.4) 



and various transition times. The curves show the fraction of fiuid A found in the 
positive extracted hemisphere, Jv. Notice that, in each case, ^ reaches a maximum 
at an intermediate transition time. These maximal fiuxes are given in Table [STT] 

Table 5.1 
Maximum percent flux ^ and corresponding r. 



Sr 


€ = V8 


^ = 7t/4 


r 


<!>(%) 


T 


<!> (%) 


{x = 0} 


3.25 


30.22 


2.7 


32.97 


{y = o} 


3.5 


64.81 


2.6 


75.68 



As r increases, the intersection curves Xq become increasingly clustered near 
the boundary of the parameter space, as can be seen in Fig. |5.5[ Eventually, this 
clustering becomes so pronounced that the curves cannot be distinguished numerically. 
This seems to imply that these regions near the boundary contribute little to the fiux; 
however, the transition map expands the small initial differences between these curves 
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Sr = {x = 0} 







- C = 


n/4 





- C = 


n/» 


o 


M.C. 


7T/4 


□ 


M.C. 


n/8 




Transition Time 



Transition Time 



Fig. 5.6. Percent composition of fluid A in the positive extracted hemisphere. Monte Carlo 
simulation results with N = 10^ are shown as the circle and square markers (see legend). The 
maximum computed flux is denoted with an * . 



by several orders of magnitude. The result is a significant contribution to the overall 
flux from the unresolved portions of the intersections. Thus, the computation of ^ for 
larger r is not numerically feasible using the action-flux formulas. However, the use 
of both forward and backward integration over the transition interval can ameliorate 
resolution problems, allowing somewhat larger r values to be reached, see Appendix [B] 
To validate our results, we again employ the Monte Carlo technique outlined in 



^4.3 Its implementation is simple since we have analytical formulas for the boundaries 
of the injection and extraction hemispheres Vq and Jv. Since Vo is simply half a ball 
of radius 1, the overall flux is estimated by 



^ 



Z ly in 



The results for A^ = 10^ are indicated in Fig. 5^ by the open circles and squares. We 
again find that the difference between the Monte Carlo and action-flux computations 
of $ is typically less than the estimated Monte Carlo error (4.15). Only for large r 



does the difference between the two become significant. Note also that the Monte 
Carlo computations are inefficient when the flux is small; most of the computational 
effort is wasted in this case since most sample trajectories do not contribute to A^^^- 



As a reflection of this, the error also grows, as indicated in (4.15). By contrast, the 



action-flux formulas remain accurate when the flux is small. For example, Monte 
Carlo simulations do not give accurate results for r < 0.5 when Sr = {x = 0} in 
Fig.m 

An advantage of the Monte Carlo method is that it appears to give reasonably 
accurate results for larger transition times than the action-flux computations. How- 
ever, this method has several limitations. If the past and future invariant regions 
had more complex boundaries than the hemispheres in our model, then initializa- 
tion of the trajectories in Vo^ and determination of whether they are advected to Jv 
would require a high-resolution representation of the lobe boundary surfaces. This 
is computationally expensive due to the exponential stretching that occurs over the 
transition interval. The same problem occurs, even when the regions have analytically 
simple definitions, if there is more than one lobe and the computation of individual 
lobe volumes is of interest. On the other hand, calculation of individual lobe volumes 
using the action- flux formulas involves no additional effort: each is computed by the 
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action-flux formulas for the particular intersection curves on its boundary. 

6. Conclusions. We have given the first quantitative analysis of transport be- 
tween Lagrangian coherent structures in aperiodic, three-dimensional flows, and have 
validated our results using Monte Carlo methods. The coherent structures for the 
transitory case are simply past- and future-invariant regions, and transport corre- 
sponds to the flux from the former to the latter. We rely on the action-flux formulas 
of ^ to provide a general framework for computing lobe volumes in n-dimensional, 
globally Liouville flows. An advantage of these formulas is that lobe volume computa- 
tions require relatively little Lagrangian information: only the orbits of codimension- 
two intersections of lobe boundary components need to be computed. Indeed, we 
found that high-resolution representations of these intersections can be obtained from 
a root-finding and continuation method, which is seeded using only a very coarse 
representation of the lobe boundary itself. An unusual aspect of the droplet model 
of ^ is that that lobes are not bounded by hyperbolic manifolds of past- or future- 
hyperbolic orbits — the bounding surfaces are simply invariant under the stationary 
vector fields P and F. 

Mixing in droplet models similar to that of ^ has been treated in [46 and, 
in the steady case, in [48]. Even though these studies computed mixing within a 
droplet, they did not address finite-time transport between coherent structures. Our 
model of smooth transitions in a serpentine mixer is perhaps more realistic than the 
instantaneous transitions used by [46], and we plan to compare our present results 
with direct simulations of the Navier-Stokes equations. 

It was observed in ^ that there is an optimal transition time that maximizes 
intradroplet transport. We would like to further investigate the effects of microchannel 
shape and choice of the injection and extraction planes on transport. Ultimately, we 
hope to find optimal channel shapes to aid in the design of efficient microfluidic mixers. 

To make this study more relevant to fluid flows, it will also be valuable to go 
beyond transport and study diffusive mixing. The inclusion of diffusive and reac- 
tive processes within the transitory framework could help in the development of a 
quantitative comparison between transport and various mixing norms. Such an im- 
proved understanding will benefit applications ranging from the design of more effi- 
cient industrial and microfluidic mixing devices to the effective large-scale recovery of 
contaminants in the ocean and atmosphere. 

Appendix A. Some Notation. Here we set out our notation, which follows, 
e.g., pLj. We denote the set of /c- forms on a manifold M by A^(M). If a G A^(M) 
and 1^1, V2, ... 14 are vector fields, then the pushfoward, R^^ of a hy R is 

(i?,a)x(yi,l^2,...,Vfe)=a«-i(,)((Di?(x))-iyi(x),...,(Di?(x))-iyfc(x)). (A.l) 

The pushforward can be applied to a vector field V as well: 

(i?*F)(x) = {DR{R-\^))V{R-\^)). (A.2) 

The pullback operator is 

R* = {R-^),. 

The interior product, or contraction, of a with V is the {k — l)-form 

zya = a(y, •,...,•). 
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Suppose that (pti,to : M x R^ ^ M is the (C^) flow of a vector field V(x, t), so that 
(ftQ^toi'^) = X, and ■^(ft^toi'^) = V{(pt^to{^)j'f^)' Then the Lie derivative with respect 
to V is the differential operator 



-^<-) -- 1 



<*o(-)- (A.3) 



The key identity for the derivative is Cartan's homotopy formula: 

jCya = iv{da) + d{ivOi). (A. 4) 

Note that C behaves "naturally" with respect to the pushfoward: 

R^Cya = Cr^vR*ol- (A. 5) 

Appendix B. Remarks on Computing Lobe Volumes. 

Here we comment on several technical aspects of our numerical implementation 
of the action-fiux formulas. We describe an efficient method for computing the inter- 
section curves X^-, discuss examples in which the action-ffux formulas may be applied 
even if the surface areas of lobe boundary components do not vanish as t ^ ±oo, 
and remark on several ways to accelerate the convergence of the integrals in (|3.5|) and 



(3.6). 



Computing the curves X^ = dVr^dJ^ri as discussed in the examples of Qand ^ 
essentially amounts to a root-finding and continuation problem. The X^ are defined 
to be zero- level-sets of some function on a 2D parameterization of dV^. This param- 
eterization is sampled with a coarse grid, neighboring grid-points that bracket the 
zero-crossing are identified, and these brackets are refined along grid lines to produce 
"seed" points that lie on 2"^. A circle of radius 5 (a pre-specified maximum Euclidean 
distance between neighboring curve points) is centered at each seed, and a ID root 
finder on the angle around the circle is used to find a new point on X^. The existence 
of this new point is guaranteed by the topology of 2"^-, provided the continued curve 
does not intersect the boundary of M. A new (5-circle is then centered at this new 
point, and the continuation is repeated until either the curve closes on itself, the angle 
between consecutive estimates of the curve tangent grows too large (this occurs when 
the curvature of X^ is large) , or the curve intersects the domain boundary. In regions 
where the curvature of X^ is large, refinement is easily performed by reducing the 
radius 5. Finally, care must be taken to ensure that a given seed does not lie on a 
previously tracked curve. 

Even if a well-defined lobe exists, and the curves X^ are computed as described 
above, it may be the case that the a-surface area of some component of the lobe 
boundary does not vanish in either direction of time under if^ and Lp^. For example 
El Uq wraps entirely around Wq{p)^ an d so i ts su rface area never vanishes under 



(/9 ; consequently, the action-fiux formulas (3.5)-(3.6) can not be used to compute the 



second integral in (4.13). We can resolve this problem by dividing Uq into subsurfaces 
whose Qf-surface areas do eventually vanish under Lp^ . For the case shown in Fig. B.l 



there are four such subsurfaces, and it is easy to see t hat U^ an d U^ collapse to p as 
t ^ — oo, and Uq and Uq collapse as t -^ +oo. Using (3.5) and (3.6), as appropriate, 
then gives 

[ a= [ If x]ds-rl[ xAds. (B.l) 

Juo J-oo \JdUl+dU^ J Jo \Jip^Q{dU^+dU^) J 
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Here we emphasize that for the second integral the past flow ip^ is to be used to evolve 
the boundaries, and the integrand is A^, the Lagrangian form for the past vector field 
P. 





Fig. B.l. (a) Subsurfaces ofUo, for the case shown in Fig.^^a), that collapse under the past 
flow (^^ as t ^' oo (blue) or t ^^ — oo (red). Arrows on Xq denote orientation with respect to a 
right-handed outward normal, (b) Subsurfaces Uq C M. Arrows denote direction of flow under ip^ . 
Blue (red) regions collapse to p in forward (backward) time. 



It is also possible to modify the application of (3.5) and (3.6) by choosing to 



evolve trajectories under a fiow that does not coincide with the transitory fiow (p. 
Indeed, as we remarked in ^ once a slice of codimension-one manifold is specified, 
Th. |3.1| is vali d for its evolution under any globally Liouville fiow. For example, the 



areas of S^ in (4.12) vanish under (p^ both as t ^ oo and as t ^ — oo. Thus, since F 
itself is a globally Liouville vector field, (3.5) may alternatively be applied with the 



replacements cp -^ cp^ and A ^ A^ to achieve the same answer as in ( 



4.12); that is. 






X^ ds. 



r(n) 



(B.2) 



We found it most efficient to choose either (4.12) or (B.2) depending upon whether 
the intersection curve is closer to /^ or /^ at t = r. In practice the intersection curves 
were closer to /^, and so we used the past integral, (B.2), for S\ and the future, 
(|4l2t,for5^. 

In addition to an appropriate choice of fiow, the convergence of the time integrals 



in (4.12) and (B.2) can be further accelerated by extrapolation. Since the intersection 



curves lie on the stable manifolds of the periodic orbits /^, they converge exponentially 
to these orbits in forward time; the contour integrals necessarily converge exponen- 
tially to zero in this same limit. Thus, the convergence of the time integral can be 
accelerated by estimating the exponential tail of its integrand (a time-dependent con- 
tour integral) using the local expansion and contraction rates, a (for the ABC case 
G = \J AC\ about the hyperbolic periodic orbits /^ (4.7). For example, the integral 
in (4.12[) can be truncated at time t to give the estimate 



f a^- f ( f x]ds-- f 

J Si Jr \Jxi J ^ Jxl 



A. 



(B.3) 



In practice we increase t until this estimate converges to some desired tolerance. 
Similar extrapolation is used to accelerate the integrals (5.15) in §51 
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Additional simplifications of the flux computation arise if the past or future vector 
fields are simple enough that explicit solutions can be found. For example, for the 
transitory ABC vector field, both P and F are integrable, and the analytical solutions 
on the separatrices are known [50 . Using these analytical formulas in the action- flux 
computations greatly reduces computation time. We use a similar simplification in 



^5.3 for the flow on the droplet boundary. When this simplification applies, the only 

numerical advection that that must be performed is over the transition interval [0,r]. 

For the computations of §5.3[ the separation of nearby trajectories affects the 



evaluation of (5.13) and (5.15) even when the intersection curves Xq are numerically 



distinguishable. Indeed the distance between neighboring points on the numerically 
computed intersection curves grows nonuniformly. Consequently, for large r, the 
resolution of Ir may be quite poor even if Xq is well-resolved. One way to ameliorate 
this effect is to use a second numerical representation of the intersection curves at 
time t = r hy solving a similar problem to (5.11) for a parameterized representation 



of <Sr. These curves can be integrated backward. Of course, in this case the resolution 
will degrade as t decreases. If we use the two representations over the first and second 



halves of the transition interval, then the accuracy of the integrals (5.13) and (5.15) 
is improved. 

Finally, if the transitory flow has symmetries, then these can be exploited to 



simplify the computations. For example, the vector field (5.7) for the microdroplet 
example, with rotations (5.10), is reversible. That is, there is a reversor Q : M ^ M 
such that 



e*F(x,t) = -F(x,r-t). 



(B.4) 



This implies that 9T = T~^6, and moreover, when Sr = %, that X^ = OXq [33 . 
Thus if we choose the numerical representations at and r to respect this symmetry, 
the forward and backward iterations over the half-transition intervals are identical. 
This is much more efficient than solving the root-finding problem (5.11 ) at both t = 
and t = r. 
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